Dynamics of delay-coupled excitable neural systems 
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We study the nonlinear dynamics of two delay-coupled neural systems each modelled by excitable 
dynamics of FitzHugh-Nagumo type and demonstrate that bistability between the stable fixed point 
and limit cycle oscillations occurs for sufficiently large delay times r and coupling strength C. As 
the mechanism for these delay-induced oscillations we identify a saddle-node bifurcation of limit 
cycles. 



1. INTRODUCTION 

The brain may be conceived as a dynamic network of 
coupled neurons [H, @, These neurons are excitable 
units which can emit spikes or bursts of electrical signals. 
In order to describe the complicated interaction between 
billions of neurons in large neural networks, the neurons 
are often lumped into highly connected sub-networks or 
synchronized sub-ensembles. Such neural populations are 
usually spatially localized and contain both excitatory 
and inhibitory neurons 

The simplest model to display features of neural inter- 
action consists of two coupled neural systems. Starting 
from this simplest network motif, larger networks can 
be built, and their effects may be studied. For exam- 
ple, starting from two interconnected reticular thalamic 
neurons with oscillatory behavior, it was shown in Q 
how more complex dynamics emerges in ring networks 
with nearest neighbors and fully reciprocal connectivity, 
or in networks organized in a two-dimensional array with 
proximal connectivity and "dense proximal" coupling in 
which every neuron connects to all other neurons within 
some radius. In another example [fj, a neural population 
was itself modeled as a small sub-network of excitable 
elements, to study hierarchically clustered organization 
of excitable elements in a network of networks. 

Most studies represent a population as an effective os- 
cillatory element that is then coupled with other pop- 
ulations to study synchronization and network effects 
@, B HI- In these studies, the basis for the emergence 
of complex network dynamics is the oscillatory behavior 
of neural elements, in spite of the fact that individual 
neural systems are usually in a stable steady state ex- 
hibiting excitable dynamics when perturbed. A reason is 
that self-sustained oscillatory behavior in the individual 
element is required to define synchronization of subsys- 
tems [l(J. One way to resolve this is to add weak Gaus- 
sian white noise to each individual element to generate 
sparse Poisson-like irregular spiking patterns, as seen in 
real neurons @, [ll], [l2j • 

In this report, we introduce a time-delayed coupling 
into the model of two excitable populations and demon- 
strate that the coupling-delay can induce sustained os- 
cillations between the two subsystems. We find a regime 
of bistability between the stable fixed point and limit 



cycle oscillations for sufficiently large delay times r and 
coupling strength C . As the mechanism for these delay- 
induced oscillations we identify a saddle-node bifurcation 
of limit cycles. 

This result suggests the use of a compound system of 
two time-delayed coupled excitable elements as a minimal 
network motif to investigate oscillatory behavior in more 
complex networks. 



2. MODEL 

We examine the delayed linear symmetric coupling of 
two identical neural populations. Each population is rep- 
resented by a simplified FitzHugh-Nagumo (FHN) sys- 
tem [l3|, [Hj], which is widely used as a paradigmatic 



model of excitable systems [15]. 

The dynamical equations are given by: 
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where the two subsystems (x\,yi) and (£2,2/2) corre- 
spond to two neuron populations. The value of the ex- 
citability parameter a determines whether the subsystem 
is excitable (a > 1) or exhibits self-sustained periodic fir- 
ing (a < 1). This is because at a = 1 the uncoupled 
systems exhibit a Hopf bifurcation, and the fixed point 
becomes an unstable focus for a < 1. e is the timescale pa- 
rameter that, if chosen to be much smaller than unity, re- 
sults in fast activator variables x\, x 2 , and slow inhibitor 
variables 2/1,2/2- Unless otherwise noted, we shall choose 
e = 0.01 and a = 1.3 for numerical simulations and thus 
restrict our analysis to parameter values where each of 
the two subsystems exhibits excitability with a stable 
fixed point. 

The interaction between the two neural populations is 
modelled as diffusive, i.e., the coupling vanishes if the 
variables identical. The coupling strength C is 

taken to be symmetric for simplicity. Further, we assume 
a finite signal propagation speed between distant neural 
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populations. This is incorporated using the delay time r. 
Note that the coupling term has the form of a classical 
diffusion term, but has lost its diffusive character through 
the introduction of a propagation delay r. More general 
delayed couplings are studied in [161 ]. 

Before we analyse Eq. (fT]) , let us briefly give a biophys- 
ical interpretation of the hitherto abstract nature of the 
dynamical variables Xi and yi. This is needed to relate 
our results to other work on population dynamics in neu- 
ral networks, in particular to studies dealing with spiking 
rates and time delays Il7l ll8l . US . 120] . These studies con- 
sider neural fields [2l|,l22|], that is, they extend over one 
or more spatial dimensions. While neural field models 
adopt the continuum limit of a network, we consider the 
minimal network motif, i. e., two discrete subsystems, but 
in both cases the biophysical reason to include delay is 
the same. 

The generic biophysical interpretation of the FHN 
model is based on a single point-like neuron and was orig- 
inally not derived from features of a neural population. 
The variable x models fast changes of the electrical po- 
tential across the membrane (spikes), and y is related to 
the gating mechanism of membrane channels [HI, [3] ■ On 
the contrary, a neural population is generically described 
by a time coarse-grained mean-field model, in which the 
dynamical variables represent averaged spiking rates [H • 
Notwithstanding, spiking rates in neural populations can 
exhibit both steady state values and relaxation oscilla- 
tions. Both states are described by the FHN mechanism 
with a > 1 and a < 1, respectively, which justifies the 
FHN model for "point-like" populations. To avoid confu- 
sion about spikes vs. spiking rates, a model of uncoupled 
subsystems as defined by Eq. (|TJ) with C — is some- 
times referred to as the Bonhoeffer-van der Pol model 
[H [Mill, for example in Ref. 



3. LINEAR STABILITY ANALYSIS 

The unique fixed point of the system is symmetric and 
is given by x* = (acj, y\, x\, y%), where x* = -a, y* = 
a 3 /3-a. 

Denoting for convenience x(t— r) = x, we can re- write 
system Eq. (fT]) as: 
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This system can be linearized around the fixed point x* 
by setting x(t) = x* + <5x(t): 
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where £ = 1 — a 2 — C. The ansatz 
$x(t) = e At u 
where u is an eigenvector of Jj? implies 
6St 
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This leads to the characteristic equation 

(l-£A + eA 2 ) 2 -(ACe- AT ) 2 =0, (7) 

which can be factorized giving 

1 - £A + eA 2 = ±ACe~ AT (8) 

This transcendental equation has infinitely many com- 
plex solutions A. Fig. [TJ shows the real parts of A for 
various values of C. As can be seen in Fig. [JJ the real 
parts of all eigenvalues are negative throughout, i. e., the 
fixed point of the coupled system remains stable for all C. 
This can be shown analytically for a > 1 by demonstrat- 
ing that no delay-induced Hopf bifurcation can occur. 
Substituting the ansatz A = ioj into Eq. ((8|) and separat- 
ing into real and imaginary parts yields for the imaginary 
part 



£ = zszCcos(ut) 



(9) 



with the Jacobian matrices J| and J* The explicit form 



This equation has no solution for a > 1 since ]£] = a 2 — 
1 + C > C, which proves that a Hopf bifurcation cannot 
occur. 



4. DELAY-INDUCED OSCILLATIONS 

The cooperative dynamics of delay-induced oscillations 
in excitable systems is inherently different from those 
of noise-induced oscillations. The introduction of noise 
terms induces sustained oscillations in each individual 
subsystem by continuously kicking these subsystems out 
of their respective rest states. Coupling then produces 
synchronisation effects between these individual oscilla- 
tors [ll], [l2]. If time-delayed feedback control (2?| is 
applied locally to one of the subsystems, the stochas- 
tic synchronisation can be tuned by varying the time- 
delay. This is in line with other work where it was 
demonstrated that such time-delayed feedback can be 
used to control the coherence and the timescales of noise- 
induced oscillations in a single FitzHugh-Nagumo system 

m mil p. 

For delayed coupling the case is entirely different. Here 
the sustained oscillations are an effect of the cooperative 
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FIG. 1: Real parts Re(A) of the eigenvalues of the fixed point 
vs. time delay r for e — 0.01, a = 1.3, and (a) C=0.1, (b) 
C=0.4, (c) C=0.8, (d) C=l, (e) C=2, (f) C=4. 



dynamics. They are generated by the delayed interaction 
between two non-oscillating stable units, and are thus 
an emergent phenomenon of the compound system. The 
bifurcation parameters for delay-induced bifurcations are 
the coupling parameters C and r. 

For large coupling delay r the oscillation is readily un- 
derstood as the two units firing alternately, each spike ini- 
tiated by the delayed signal of the remote system. Since 
the signal of one system is transmitted to the other and 
then back, the oscillation in one system must have a pe- 
riod T of approximately 2r, and be phase shifted by T/2 
with respect to the other system. This is visible in the 
time series (Fig. (2fa),(b)) and in the phase portraits of 
the activators (Fig. Hfc)) and inhibitors (Fig. [5] (d)). 

Since the two subsystems have identical but T/2- 
shifted profiles Xi{t) and yi(t), the coupling terms 
C[xi(t— t)— Xj(t)] (i ^ j) would vanish for all t if the 
oscillation period T were precisely 2t. Then the periodic 
orbit (Fig. [3^b)) would also be a solution for the uncou- 
pled subsystems (C = 0). This cannot be the case for 
the excitable regime of the FHN systems (a>l) because 
periodic orbits do not exist there. Hence it follows that 
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FIG. 2: Delay-induced oscillations, (a), (b): Time series of 
both subsystems (red solid lines: activator Xi, green solid 
lines: inhibitor yi\ black dashed lines: fixed point values of 
activator and inhibitor), (c), (d): Phase portraits of activa- 
tors (c) and inhibitors (d). Parameters: e = 0.01, a — 1.3, 
C = 0.5, t = 3. 



the oscillation period T must be of the form T = 2(t + S) , 
where 5 ^ is the effective time shift responsible for the 
non-vanishing coupling term. The phase portraits of the 
activators Xi in the planes (xi(t), x^it — r)) deviate from 
the bisector (Fig. SJa)) due to the fact that Sy^O. 

To obtain a clear picture of the timescales involved 
in the dynamics, we have computed the excursion times 
along the segments of the phase-space trajectory. The 
start and end points of the different segments of the 
trajectory (colored dots A to D) are defined to corre- 
spond to the time steps when the trajectory has left or 
entered a neighborhood Axi (here Axi = 0.01) of the 
stable branches of the xi-nullclinc. 

In the example shown in Fig. [3] (e = 0.01, a —1.3, C = 
0.5, r = 3) the value of the effective time shift is <5 — 0.012. 
This value is about one third of the fast transition times 
between the stable branches of the xi-nullcline, Ati = 
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FIG. 3: Blow-up of delay-induced oscillation in the first sub- 
system: (a) time series and (b) phase portrait (xi,yi). The 
four different stages of the limit cycle are separated by col- 
ored dots A, B, C, D. Parameters: e = 0.01, a — 1.3, C = 0.5, 
t = 3. 
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0.041 and Ai 3 = 0.038, respectively. Aii is the rise-time 
of the spike, i.e., the time that elapses between leaving 
the fixed point (black dot A) and crossing the right stable 
branch of the ^i-nullcline at B (green dot, Fig. (2Kb)). 
Ai 3 is the drop-time, i.e., the duration of the jump back 
from the right to the left stable branch of the xi-nuhclinc 
(between blue and magenta dots, C — > D). The slow 
parts of the trajectory, occurring on the right (B — > C) 
and left (D — > A) stable branches of the xi-nullcline, have 
a duration Ai 2 = 0.357 and At 4 = 5.588, respectively. The 
total oscillation period is thus T= 6.024. 

In Fig. O the phase portraits of the first of the two 
delay-coupled subsystems are shown for different ex- 
citability parameters a and delay times r. The top panel 
(a,b) corresponds to excitability parameters far from the 
Hopf bifurcation of the uncoupled system, which occurs 
at a = l. The bottom panel (c,d) corresponds to values 
of a close to the Hopf bifurcation. 

In the case of a = 1.3 and r = 3 (Fig. EJa)), the oscil- 
lation period T — 6.024 is large enough for the two sub- 
systems to nearly approach the fixed point before being 
perturbed again by the remote signal. Note that Fig. 0(a) 
corresponds to Fig.[3£b). If the delay time becomes much 
smaller, e.g., for r = 0.8 (Fig. 0(b)), the excitatory spike 
of the other subsystem arrives while the first system is 
still in the refractory phase, so that it cannot complete 
the return D — > A to the fixed point. The times spent 
in the different stages of the limit cycle for r = 0.8 are 
Aii = 0.087, Ai 2 = 0.111, Ai 3 = 0.085, and Ai 4 = 1.355, 
and the oscillation period is T= 1.637. The effective time 
shift 5 = 0.018 for r = 0.8 is larger than for r = 3. 

We now investigate the transition from large to small 
delay times r close to the Hopf bifurcation at a = 1.05. 
Again we choose r = 3 (Fig. [SJc)) and find T = 6.018 
and 8 = 0.009. The times spent in the different phases 
of the oscillation period T are Aii = 0.043, Ai 2 = 0.329, 
Ai 3 = 0.186, and Ai 4 = 5.461. For the shorter delay 
time t = 0.8 (Fig. E^d)) we find T = 1.630, S = 0.015, 
Aii = 0.072, Ai 2 = 0.240, Ai 3 = 0.070, and Ai 4 = 1.249. 
We find the same pattern as far from the Hopf bifurca- 
tion. The effective time shift S is smaller if r becomes 
larger and the system can approach the fixed point. Fur- 
thermore, if the system is close to the Hopf bifurcation 
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FIG. 4: (a) Phase portrait of delay-coupled excitable sys- 
tem in the plane (xi(t),X2(t — r)) with positions on the orbit 
marked as colored dots A, B, C, D as in Fig. [3] (b) Time series 
of (x2{t — t) — x\(t)). The inset shows a longer time series. 
Parameters: e = 0.01, a — 1.3, C = 0.5, r = 3. 
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FIG. 5: Phase portraits of delay-coupled excitable system 
(xi,yi) for different excitability parameters a and delay times 
t (trajectories: solid blue, nullclines: dashed black), (a) 
a = 1.3, t = 3, (b) a = 1.3, r = 0.8, (c) a = 1.05, r = 3, (d) 
a = 1.05, r = 0.8. Other parameters: e = 0.01, C = 0.5. 



and t is sufficiently large so that the fixed point is very 
closely approached, S becomes much smaller than far 
from the Hopf bifurcation, and tends to zero for a — * 1. 
The reason is that when the excitatory spike arrives in 
the first subsystem, there is a turn-on delay 5 before the 
first subsystem emits a spike. This is because the trajec- 
tory has to cross the middle branch of the xi-nullclinc 
due to this upstream impulse. This section of the tra- 
jectory becomes the smaller, the closer the fixed point A 
(at x\ = —a) is to the minimum of the x\ nullcline (at 
x\ = — 1). This also explains the origin of the time-shift 
S. 

Finally, we shall investigate the question whether the 
system exhibits bistability between the fixed point and 
the limit cycle oscillation for all values of r. In Fig. [B]the 
regime of oscillations is shown in the parameter plane 
of the coupling strength C and coupling delay r. The 
oscillation period is color coded. The boundary of this 
colored region is given by the minimum coupling delay 
Tmin as a. function of C. For large coupling strength r m , n 
is almost independent of C; with decreasing C it sharply 
increases, and at some minimum C no oscillations exist 
at all. At the boundary, the oscillation sets in with fi- 
nite frequency and amplitude as can be seen in the inset 
of Fig. [6] which shows a cut of the parameter plane at 
C = 0.8. The oscillation period increases linearly with 
t. The mechanism that generates the oscillation is a 
saddle-node bifurcation of limit cycles (see Fig. [5] inset 
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FIG. 6: Regime of oscillations in the (r, C) parameter plane 
for initial conditions corresponding to single-pulse-excitation 
in one system. The oscillation period T is color coded. The 
transition between black and color marks the bifurcation line. 
Inset (a) shows the oscillation period vs. time delay r in a cut 
at C — 0.8. Parameters: e = 0.01, a = 1.3. Inset (b): schematic 
plot of the saddle-node bifurcation of a stable (red solid line) 
and unstable (blue dashed) limit cycle. The maximal oscilla- 
tion amplitude is plotted vs. the delay time r and the stable 
fixed point is plotted as a solid black line. The gray back- 
ground marks the bistable region. 



(b)), creating a pair of a stable and an unstable limit 
cycle. The unstable limit cycle can be visualized from 
transients starting at appropriate initial conditions (see 



CONCLUSION 



We have shown that delayed coupling can induce peri- 
odic spiking in a compound system of two coupled neuron 
populations if the delay and the coupling strength are suf- 
ficiently large. Bistability of a fixed point and limit cycle 
oscillations occur even though the single excitable ele- 
ment displays only a stable fixed point. The two neural 
populations oscillate with a phase lag of it. 
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FIG. 7: Unstable delay-induced limit cycle, (a), (b): Time 
series of both subsystems (red solid lines: activator Xi, green 
solid lines: inhibitor yi\ black dashed lines: fixed point value 
of activator and inhibitor). (c),(d): phase portraits of (x\,y\) 
(top) and (22,3/2) (bottom). The unstable small-amplitude 
limit cycle is visualized by the initial part of the transient 
approaching the stable fixed point. Parameters: e = 0.01, 
a = 1.3, C = 0.5, r = 3. 
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